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I Abstract 

u 

We  study  some  of  the  numerical  properties  of  the  nested 
decomposition  algorithm  of  Ho  and  Manne.  In  particular  we  seek 
to  show  how  well  developed  theory  in  the  area  of  computational 
linear  algebra,  due  primarily  to  J.H.  Wilkinson,  carries  over  to 
linear  programming  and  yields  useful  insight  into  the  behavior 
of  algorithms  in  this  area. 
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ERROR  PROPAGATION  AND  SOLUTION  RECONSTRUCTION  IN  NESTED  DECOMPOSITION 

by 

L.  Nazareth 


1.  Introduction 

The  Idea  of  a nested  decomposition  algorithm  dates  back  to  the 
fundamental  paper  of  Dantzig  and  Wolfe  [I960].  Many  authors  have  since 
contributed  to  its  development  including  Dantzig  [1963,  Ch.  23],  Cobb  and 
Cord  [1967],  Glassey  [1973],  Ho  and  Manne  [1974],  Ho  [1974a,  1974b], 

Kallio  [1975].  The  major  credit  for  the  development  of  this  algorithm 
goes  to  Ho  [1974a, b],  who  put  it  on  a sound  algorithmic  footing,  implemented 
the  algorithm  and  showed  it  to  be  a workable  technique  for  solving  LP 
problems  with  staircase  structire,  of  the  form: 


minimize  £ 
t-1 


c.  x 


subject  to  A^x^ 


4 


(l.D 


Bt-A-i  + At*t " it 


x.  > 0 


for  t * 1, . . . ,T 


where  x is  n„  x l A_  is  m x n , 13  ■ in  , , x n , c is  1 x n , and 

— t t t t t t t*Tl  t — t t 

d^  is  mt  x 1 in  dimension.  We  shall  call  the  above  problem  LP. 

— t c 

In  this  paper  we  look  at  Ho's  algorithm  from  a numerical  standpoint. 
This  paper  is  organized  into  five  sections  as  follows: 
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In  Section  2 we  provide  background  by  briefly  describing  the 
algorithm  and  relating  practical  experience  obtained  by  M.  Aganagic*  [1977] 
who  extended  the  code  of  Ho  and  ran  it  on  a version  of  the  PILOT  model  of 
Dantzig  and  Parikh  [1975].  It  is  this  experimentation  which  motivated  the 
research  described  here. 

In  Section  3 we  give  examples  which  illustrate  numerical  difficulties 
which  can  arise,  and  we  discuss  certain  numerical  properties  of  this  algorithm. 

In  Section  4 we  use  some  of  the  results  of  backward  error  analysis 
and  perturbation  theory  of  computational  linear  algebra,  Wilkinson  [1965], 
in  order  to  study  the  propagation  of  error  from  stage  to  stage. 

In  Section  5 we  examine  further  the  reconstruction  of  the  solution 
in  nested  decomposition,  in  the  presence  of  numerical  error.  We  propose 
an  alternative  method  for  doing  this. 

Finally  in  Section  6 we  make  some  recommendations  about  the  imple- 
mentation of  nested  decomposition  and  describe  some  lessons  we  have  learned 
from  our  research,  about  the  development  of  structured  LP  algorithms  and 
codes. 

Our  main  aim  in  this  paper  is  to  look  at  LP  algorithms  based  upon 
the  decomposition  principle,  from  a numerical  standpoint  and  to  show  that 
such  algorithms  have  some  very  Interesting  numerical  properties.  A complete 
investigation  is,  of  course,  well  beyond  the  scope  of  this  paper. 


Ph.D.  student.  Operations  Research  Department,  Stanford  University,  California. 
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2.  Background  I ] 

< 

2.1.  Brief  description  of  algorithm 

We  assume  familiarity  with  the  nested  decomposition  method  for  solving 
systems  of  the  form  (1.1)  as  described  In  Ho  [1974a].  We  summarize  this 
by  the  schemata  shown  in  Fig.  1,  for  a four-stage  problem. 

. 

Notes  on  Figure  1. 

(a)  We  distinguish  extreme  points  and  corresponding  proposals  by  super- 
scripts, to  emphasize  that  we  are  dealing  with  vectors. 

1 

(b)  The  matrix  of  generated  proposals  is  denoted  by  . In  order  to 
emphasize  the  distinction  between  a master  problem  and  the  correspond- 

j i 

ing  restricted  master  problem  we  append  the  symbol 

i I Ml  1 

(c)  The  convexity  row  is  denoted  by 

(d)  The  objective  row  for  a subproblem  comes  from  the  pricing  out  of 
proposals  generated  by  the  subproblem,  as  explained  in  the  right-hand 
column. 

2.2.  Reconstruction  of  the  solution 

This  has  again  been  described  in  detail  in  Ho  [1974a].  He  gives 
two  omthods.  In  the  first  method  (Method  I),  all  extreme  points  developed 
for  each  subproblem  are  kept.  Let  (X^,  x^)  be  the  optimal  solution  of 
S/P  4.  In  Fig.  1, 


In  the  second  method  (Method  II) , it  is  unnecessary  to  keep  all 
extreme  points.  Instead  they  can  be  regenerated,  if  necessary,  in  the  follow- 
ing reconstruction  algorithm.  Let  us  again  consider  the  four-stage  problem. 

* 

Fix  x^  at  x^  and  treat  the  problem  as  a three-stage  decomposition  in 
variables  x^,  £2*  and  x^.  Following  Figure  1,  the  last  stage  would 
then  be  of  the  form 


I / / 


-3 


>S/P  3 


“ V4  ' 


} Variables 


and  the  first  and  second  stage  would  be  as  in  Fig.  1.  When  this  problem  is 
solved  x-j  will  be  known,  and  let  us  say  its  optimal  value  is  x^.  We 
can  now  fix  x^  and  x^  at  their  optimal  values,  i.e.  the  last  set  of 
equations  effectively  drops  out.  We  now  solve  the  first  three  sets  of 
equations  as  a two-stage  decomposition  in  variables  jx^  and  x^*  This 
process  is  continued  until  the  optimal  value  of  all  variables  is  known. 

Note:  Reconstructed  solutions  are  not,  in  general,  basic  solutions 
for  the  original  LP. 
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2.2.  Practical  experience 


By  extending  the  code  of  Ho  and  running  it  on  a version  of  the  PILOT 
Model,  of  Dantzig  and  Parikh  [1975],  M.  Aganagic  [1977]  obtained  a great  deal 
of  valuable  experience  with  the  algorithm.  He  found 

• that  master/subproblems  could  become  badly  scaled 

• that  in  certain  cases,  propagation  of  error  from  stage  to  stage  was 
substantial 

• that  the  technique  for  reconstructing  the  solution  could  fail. 

However  the  code  used  was  an  experimental  one,  and  the  LP  problem 

itself  had  coefficients  which  differed  widely  in  magnitude.  It  is  therefore 
unclear  whether  these  difficulties  are  inherent  in  the  nested  decomposition 
algorithm,  or  whether  some  of  them  are  the  result  of  numerical  instabilities 
in  the  code  (associated  with  use  of  the  product  form  of  the  inverse)  or 
inherent  in  the  LP  problem  being  solved.  Hence  the  analysis  described  here. 
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3.  Examples  and  Diacuslon 


We  now  discuss  each  of  the  difficulties  mentioned  in  Section  2.2, 
mainly  through  numerical  examples. 

3.1.  Scaling. 

An  LP  problem  which  is  well  scaled  for  the  simplex  algorithm,  can 
be  badly  scaled  for  the  decomposition  algorithm.  Consider  the  problem 

minimize  x + y 

subject  to  x + ey  > 0 
x - ey  > 0 
x - ey  <_  1 
x + ty  < 1 

y < i 
-y  < i 

where  e is  small. 

If  the  simplex  algorithm  is  applied  to  this  problem,  every  basis  is 
well-conditioned  (i.e  has  a reasonable  condition  number),  and  small  variations 
in  the  coefficients  will  produce  small  variations  in  the  extreme  points  of 
the  problem. 

Suppose  however  we  decompose  this  LP  into  a master  and  subproblem  as 
indicated  above . The  extreme  points  of  the  subproblem  are 
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CM:)- 

and  the  master  problem  becomes 

(l/2e)A3  + (-l/2e)A4  + s3  - 1 

(-l/2e)\3  + ( l/2e)A4  + s2  = 1 

A.  + A„  + A_  + A.  - 1 

12  3 4 

2 8i»  82  _>  0 

This  problem  has  basis  matrices  with  a large  condition  number 
~ 0(l/e)  and  the  structural  (non-slack)  columns  differ  widely  in  magnitudes 
so  that  their  reduced  costs  can  give  misleading  information  about  the  value 
of  introducing  a particular  column  into  the  basis.  In  this  sense,  see 
Tomlin  [1975],  the  master  is  badly  scaled. 


3.2.  Error  propagation 

Consider  the  two-stage  problem 


Vi  ■ *i  } 

Bl-1  + *2-2  “ -2  } 

*1*  *2  - 0 


S/P  rows 
Master  rows 


(3 
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Let  X}  be  the  true  extreme  points  of  the  S/P  and  xj  the 
computed  extreme  points.  If  a stable  algorithm  is  used,  then  by  the  back- 
ward error  analysis  of  Wilkinson  [1965],  each  x|  is  the  exact  solution 


of  a perturbed  system 


(Aj  + 


(3.2b) 


where  a|  is  the  basis  matrix  corresponding  to  and  is  a small 

perturbation.  For  example,  if  we  assume  that  all  elements  of  satisfy 

I (Al>iml  — 1 ’ tllat  Partia^  pivoting  in  the  LU  factorization  of  a|  is  used 
and  that  |u  | < 1,  then  ||E^||  < 3*n*2  1,  when  t,  is  the  number  of  bits 

in  a floating  point  word  (see  Wilkinson  [1965]). 

Note  that  x^  can  be  extremely  different  from  xj»  this  being 
determined  by  the  condition  number  of  a|.  We  can  show  by  perturbation 
theory,  Wilkinson  [1965],  that 


y - 31  MiiE%yii)  . . . 

— < — J 3 7—  where  k - ||aJ ||||(AJ)~1|| 

llxjll  (1  - kj  || E^  || / II A^ || ) J 1 1 

If  the  exact  solution  to  (3.2a)  is  x^  - £ X x^ , 7 X • 1,  and 

J J J J 

computed  solution  x^  « ^ ^j^l»  Ij  Vj  ■ then 

llic  ’ Sjl  *11 1 *j<Ai  + - l Xj  (aJ)'”H1|| 

assuming  (A^  + E^)  and  a|  are  invertible, 

- ||£  Xj(I  + Aj  V)"1^)"1^  - l XJ(A})~1b1H 


' ' • * 


Now,  for  any  square  matrix  X,  (I  + X)-1  - (I  - X + X2  - X3  +•••),  provided 
II X||  < 1.  If  we  assume,  therefore,  that  ||(Aj)-3E3||  < 1 


X-  - x. 


- - ^)Aj  bj  + l Xj(-(A3)"1Ej)(I-(Aj)'1Ej+---)(Aj)~1b1l| 

< III " xj)(A3)'1b1ll  + III  ^((A^rVxi  + (Aj)"1Ej)_1(A^)~1b1» 
^ J 


Assume  further  that  A j*  0 <»»>  A + 0. 

J j 


Hxc-xtll  < 


j(^  " )XW  -1"  + 11  ^XJ(Ai)~1-i)(a^)(a1)"1eJ(i+(a^ 


J-V)-1! 


11 W 


< I 


\ \ l|A^ II II ( AJ ) -1  H ( I|EJ n / II H > 

XJ  / Ijd  - II^IIII(aJ)"1||(||ej||/||aJ||) 


This  bound  is  not  an  encouraging  one,  if  ||a|||||(A^)_1|| 


differs  substantially  from 


is  large  or 


Let  us  however  now  adopt  the  viewpoint  of  backward  error  analysis. 
Ignoring,  for  the  moment,  errors  in  optimizing  the  master  (see  Section  4) 
then  we  are,  in  effect,  exactly  solving  a system  of  the  form  (3.2b),  in 
which  each  element  of  A^  has  a small  uncertainty  with  a known  bound. 


NOTE  that  since  E^  varies  with  A3,  we  cannot  have  A^  perturbed 
by  a fixed  matrix,  i.e.  by  a static  perturbation.  We  shall  use  the  term 


dynamic  perturbation  to  distinguish  our  case. 


Now,  if  we  assume  that  the  solution  of  the  LP  is  stable,  then 
we  would  expect  - x^l  to  be  small.  This  implies  a correlation  of 

error,  which  the  following  example  exhibits: 


minimize 
subject  to 


*!  + (l+e1)x3 

X1  + *2  + d+e2)x3 


x3  + x4 


x±  > 0 


a2 

a3 

a4 


S/P 


Master 


(3.2c) 


There  are  only  two  feasible  extreme  points  of  S/P,  and  these  are 


when  k ■ (a^-a^)/ (e2-c^)  and  where  we  assume  that  and  e 2 are 

small,  e2  > eL  and  a2  > a1  > a2(l  + 6^/(1  + e2>.  The  master  problem  is: 


minimize  x. 

4 


i 
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This  has  the  solution 


x4  “ a4 


^2  “ ~ a4)  /k 

A - 1 - (a3  - a^)/k 


where  we  assume  a^  > a^  X 0 and  chosen  so  that  X2  ^ 0«  The  optimal 

1 2 

solution  xQpt  is  given  by  + ^2—1’  becomes 

a1  - (1+e^  (a3-a^,) 

(82-32)  - (a3~a^)  (£2^2) 

^opt 

a3'a4 


Solving  the  system  (3.2c)  as  an  LP  also  yields  the  unique  (optimal)  solution 
and  the  associated  4x4  basis  is  well  conditioned.  This  implies  that  x0pt 
is  relatively  insensitive  to  small  changes  in  the  matrix  or  r.h.s.  elements. 

AAA  * 

Suppose  that  a2»  a2,  tj,  t2  are  obtained  by  truncating  82 , a2,  £>2  and  e2 
which  are  not  machine  representable  numbers.  Suppose  also  that  the  above 
derivation  is  carried  out  with  &2»  a2’  ^1  and  e2  used  in  place  of  ai»  a2* 
£2*  and  e2,  and  denote  the  quantities  corresponding  to  k,  X2  and  X2  by 

AAA  ^ 

k,  X2  and  X2>  Then  the  quantity  k which  determines  ^ can  be 

A A A A A 

drastically  different  from  k - / <£  2~z.y) . This  is  because  the 

2 

basis  which  determines  ls  Ul-conditioned.  The  corresponding  quantities 

Xj  and  X2  will  also  be  quite  different  from  X2  and  X2-  However  the 

1 2 

errors  in  Xj  and  X2  are  correlated  with  the  errors  in  and  Xj,  so 

that  x _ does  not  change  drastically. 

~^opt 


3.3.  Reconstruction 


Finally  we  give  an  example  which  illustrates  difficulties  associated 
with  reconstruction  of  the  solution.  We  wish  to  show  that  Method  2 of 
Section  2.2  can  result  in  numerical  instability. 

Consider  the  problem: 


maximize 


subject  to  x^  + (l-te^)Xj  - 


x2  + (1+e  2>x3 


“ b] 
- b. 


S/P 


(3.3a) 


Xi  + x2  + 2(l4e  3>Xj  + x4  - b3  j Master 

x±  > 0 


There  are  only  two  feasible  extreme  points  and  these  are  given  by 


/ b -b2(Uc ,)/(l+e2) 

X!  - ( b2 

• i'l  • 

\o  j 

\ b2/(l4e2) 

(3.3b) 


1 2 

where  we  assume  that  b^  and  b2  are  chosen  so  that  — ]_  - All 

bases  corresponding  to  these  are  well-conditioned.  The  master  problem 


is  given  by 


maximize 


subject  to 


(b2  + b2)x1  + 


b (I4e  ) 2(He  )b 
bl  ' "(He  J + (1-fcJ 


X2+X4  " b3 


(3.3c) 


Xl  + 


X2  * 1 


V 2* 


x.  > 0 
4 — 


Assume  (2e^  - e^)  < e2.  Then  the  coefficient  of  X2  < coefficient  of  Xj 
and  so  the  optimal  solution  is  given  by  X^  ■ 0,  X2  ■ 1 and  determined 

by  (3.3c).  Then  the  optimal  solution,  by  Method  I of  Section  2.2, 


x _ 1 
—opt 


b1  - b2(l-K:1)/(l4c2) 


b2/a*2) 


(b1  + b2(l  + 2e3  - El)\ 

\ Tr+-2i ) 


(The  basis  corresponding  to  variables  x^,  x^  and  x^  in  the  original  LP 

is  well-conditioned  and  the  solution  x . can  be  alternatively  obtained 
— ■■—..*■  opt 

from  this  basis.)  If  Method  II  (cf.  Section  2.2)  is  however  used,  we  must 
then  solve  the  system 
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1 0 I-K.  1 

x. 

r hi 

1 

1 

i 

0 1 1+e  2 

X2 

= 

b2 

(1  + 2(e  -e  ))b. 

_ 1 1 2 (1+e  3>_ 

- x3  - 

i_  J 1 z 

L bl  (1  + e„) 

This  is  an  ill-conditioned  problem,  and  thus  numerical  error  could  result 
in  the  reconstructed  solution  being  substantially  different  from  the  true 
we 11- conditioned  optimal  solution,  even  possibly  infeasible.  We  return 
to  the  problem  of  reconstruction  in  Section  5. 
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4. 


Analysis  of  Error  Propagation 

In  this  section  we  use  some  of  the  results  In  Wilkinson  [1965],  « 

Chap.  4 on  backward  error  analysis  and  perturbation  theory  of  linear  systems 
of  equations,  in  order  to  study  the  propagation  of  error  from  stage  to 
stage.  A complete  error  analysis  is  beyond  the  scope  of  this  paper. 


4.1.  In  backward  error  analysis  one  seeks  to  show  that  the  computed  solution 
x of  a problem,  say  P,  is  the  exact  solution  of  a problem  obtained  by 
perturbing  P.  In  the  nested  decomposition  algorithm,  a typical  restricted 
subproblem  is  of  the  form  (see  Fig.  1) 


m rows:  Q A + Ax  « b 

t L t C — t L 


convexity  rows:  e A 


- 1 


S/P  - t 


(4.1a) 


*t»  -*t  - ° 


where  Qt  - and  X^  ^ is  the  matrix  whose  columns  are  computed 

T 

extreme  points  and  rays  of  the  previous  subproblem  S/P  - (t-1).  e^  is  a 

c 

row  vector  whose  j’th  element  is  1 if  the  corresponding  colunn  of  i® 

an  extreme  point  and  0 if  it  is  an  extreme  ray.  A typical  basis  of 
S/P  - t consists  of  a set  of  mfc  + 1 columns,  and  is  of  the  form 


(4.1b) 
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m • - m M 


* A * «p  Q J 

where  At  and  ef  denote  columns  of  xt_^»  At  and  ®.t« 

Computed  extreme  points  and  rays  for  the  next  subproblem 
S/P  - (t  + 1)  are  exact  solutions  of 


(6*  + 6$)*^ 


(4.1c) 


e.g.  if  (i) 

tc^ijl  < 1 

(ii) 

partial  pivoting  in 

the  LU  factorization 

of  is  used 

(iii) 

'“u'i1 

-t 

Then  HfiOll^  1 3 

(mt  +1)2  , see 

Wilkinson  [1965]  and 

(4. Id) 


Under  what  circumstances  can  we  cast  errors  in  <8  back  into  the  data  of 
the  original  problem  and  into  the  convexity  row? 

Let  us  partition  6®  in  the  same  way  tha':  Q is  partitioned 
in  (4.1b).  Then 


£<8 


■«0u 

5®12  ' 

**21 

tm 

6 ^22 

(4. le) 


The  errors  £<^2  can  attr^ute^  to  At*  in  (4. lb) and  [6 ^21  ^ ^ ®22 ^ 
to  the  convexity  row  in  (4.1b).  In  addition  let  us  seek  £Bt_^  such  that 
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<Vl  + “t-X-l  - Vl*t-1  + s«u 


or 


5Bt-ixt-i  * 


11 


A strong  assumption  which  ensures  this  is  that  every  matrix  X 
is  a feasible  basis  has  full  column  rank  and  that 


t-1 


+ T 

X®.!  - (X®_x  X®.!)’1  ^ satisfies 


»cT 


Xt-lH  i k • 


where  Z denotes  the  generalized  inverse  of  Z.  Then 


|l«Bt_1IL  < Mn\\J  XJL^  < 3k(mt  + 1)2 


-t 


when  (4. Id)  holds. 

A c 

Note  that  6Bt_i  varies  with  Xt_^.  We  can  however  say  that  when 

the  above  assumptions  hold,  the  computed  solution  is  the  exact  solution 

obtained  by  nested  decomposition,  applied  to  an  LP  (1.1),  whose  matrices 

At  have  a dynamic  perturbation  bounded  by  3(mt  + 1)2  , and  whose  matrices 

”*■1 

Bt  have  a dynamic  perturbation  bounded  by  3* k* (mt+^  + 1)2  . Further  each 

convexity  row  in  the  nested  decomposition  has  a dynamic  perturbation  bounded 
by  3* (mt  + 1)*2 

The  above  assumptions  are  very  restrictive.  Less  restrictive 
assumptions  would  require  6$  ^ to  be  cast  into  errors  in  the  data  of 
previous  stages.  Also  in  order  to  circumvent  having  to  explicitly  introduce 
error  into  the  convexity  row,  we  could  force  it  to  be  satisfied,  by  using 
a method  similar  to  that  used  in  the  GUB  algorithm  Dan tzlg  and  Van  Slyke  [1967], 
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4.2.  Perturbation  theory  complements  backward  error  analysis,  see  Wilkinson 
[1965].  We  study  the  following  questions: 

Define  LP^  by  perturbing  each  matrix  A£  by  6At  and  each  Bt 
by  6Bt  in  (1.1),  where  6At  or  6Bt  are  small  fixed  (static)  perturbations. 
If  the  matrices  Qt  in  the  decomposition  of  LP  are  in  consequence  perturbed 
to  Q*  we  ask:  Will  each  column  of  be  close  to  the  corresponding 

columns  of  Qt? 

We  make  the  following  strong  assumptions: 

Assumption  4.2:  Every  basis  of  S/P  - t (see  (4.1a))  is  well-conditioned 

and 

ll<2tllll$;1H  < c (4.2c) 

Assumption  4.3:  The  subproblems  in  the  decomposition  of  LP^  match  those  in 
the  decomposition  of  LP  in  the  following  sense:  each  has  the  same  number 
of  columns;  furthermore,  if  a basis  & t is  feasible  for  SP  - t in  the  nested 
decomposition  of  LP,  the  basis  for  the  corresponding  set  of  columns 

in  SP^  - t,  in  the  decomposition  of  Lp\  is  also  feasible. 

Given  a basis  for  subproblem  t of  LP,  let  + 6^ 

for  some  6($t  be  the  corresponding  basis  for  the  subproblem  of  LP*.  Let 
the  two  feasible  solutions  be  x^  and  x^  + Sx^,  respectively. 

Thus 

«txt  * bt 

+ *St)(»t  + «»t>  ■ *>t 
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As  in  Section  3.2 


l|6xtl|  c||60t||/||0t|| 

n^r  1 (i  - cii6®tn/ii^tii) 

where  c is  defined  in  Assumption  4.2. 

The  proposal  transmitted  to  SP  - (t+1)  is  qt+^  * Btxt  and  to 
SP1  - (t+1)  is  q^+1  - (Bt  + 6Bt)(xt  + 6xt).  So 

HqJ+l  " qt+1l!  < H6Btxt||  + ||Bt6xt|| 


l|6Btl|||xtH  l|Btl|c(||6  0tH/||<fitl|)  llxjl 
H'BtxJ  + (1  - c(||6<rtI/1Wtll»  lfBtxt|| 

If  we  assume  that  each  qt+i(q^+^)  in  a Basis  is 

of  the  same  order  of  magnitude,  the  left  hand  side  of  (4.2c)  is  also  an 
estimate  of  \\6&  t+1ll/ll$t+1IU 

If  we  therefore  let  l|60jJ|/ 11(3^11  <.2  \ and  note  that 

llBtxtll  i llBtl| ||xt ||  so  that 


ik 


t+1 


- q 


t+i" 


"Ik 


t+i1 


c(||6fitl|/||fltH) 
a - c(\\6&t\\/wt\\)) 


underestimates  the  2nd  term  in  (4.2c),  then  it  is  readily  seen  by  induction 
chat  the  leading  term  in  the  right  hand  side  of  (4.2c)  contains  the  factor 


(c  2 ).  This  suggests  that  the  columns  of  Qt  could  depart  substantially 

from  the  columns  of  Qt  when  t is  large.  Note  however  that  the  solution 
of  LP*  will  be  very  close  to  the  solution  of  LP  if  the  optimal  basis  is 
stable.  Again,  as  in  Section  3.2,  we  have  the  effect  of  correlation  of  errors. 


5.  Reconstruction  of  the  Solution 


In  Section  2.3  we  showed  an  example  of  how  numerical  difficulties 
can  arise  with  reconstruction  of  the  solution.  In  this  section  we  discuss 
further  difficulties  associated  with  solution  reconstruction. 


5.1.  Returning  to  the  two-stage  problem  (3.2a)  of  Section  3.2,  suppose  now 
that  computed  proposals  x*,  which  are  exact  solutions  of  (3.2b),  are  such 
that  IIE-*  ||  is  no  longer  small — e.g.  where  an  unstable  algorithm  is  used 

to  solve  the  system  of  equations  A^x^  “ b^ . In  this  case,  the  optimal 

* * 1 
solution,  say  (x^,X2)  a linear  programming  problem  LP  , obtained  by 

dynamically  perturbing  (3.2a)  , may  be  such  that  a no  vector  x.^  > 0 for 

which  (x^,X£)  is  feasible  for  LP  (3.2a).  (Recall  that  the  (dynamic) 

perturbation  is  determined  by  the  pattern  of  decomposition.)  Thus  method 

II  of  Section  2.2,  which  solves  directly  the  LP 


minimize 


subject  to 


(5.1a) 


could  fail.  A similar  situation  could  arise  when  the  assumptions  of  Section  4.2 
do  not  hold,  and  it  Is  not  possible  to  account  for  computational  errors  in 
the  master  problem  in  terms  of  small  (dynamic)  perturbations  in  the  original 
data  A^,  A2,  B^.  Note  that  Method  I will  succeed,  in  the  sense  that  It 
reconstructs  a solution  which  is  feasible  for  the  perturbed  LP.  Thus  in  the 
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presence  of  numerical  error  Method  II  does  not  adequately  reproduce  Method  I. 
Similar  arguments  also  hold  for  a t stage  problem.  It  is  also  clear  from 
these  arguments  how  crucial  it  is  to  use  a stable  implementation  of  the 
simplex  method  for  solving  the  subproblems. 

The  second  point  we  would  like  to  make  is  that  loss  of  feasibility 
occurs  in  a particular  way,  as  we  now  show.  Consider  a three-stage  problem 
of  the  form  (1.1).  As  summarized  in  Fig.  1,  the  true  extreme  points  of  the 
first  stage  are  denoted  by  x^»  and  we  assume  that  the  corresponding  sub- 
problem is  bounded.  The  true  extreme  points  of  the  second  stage  also  assumed 


bounded  are 


and  we  shall  assume  that  these  are  deliberately  perturbed  to 


IV  J . The  third  stage  master  is  then  given  by 


I 


I 

J 

Let  the  optimal  solution  obtained  by  solving  (5.1b)  exactly  be 


Then  by  Method  I 


Note  now  that  (x^  x2,  x3>  satisfies  stage  1 and  stage  3 
constraints  exactly.  Only  stage  2 constraints  are  violated.  (We  have 
assumed  that  2.^  *2k  “ **  a8ain  demonstrating  the  values  of  solving  sub- 
problems by  a GUB  type  algorithm.)  Error  has  however  propogated  to  stage  3 
in  the  sense  that  allocation  of  resources  to  each  stage  is  affected  and 
hence  so  is  the  computed  optimal  solution.  However  feasibility  of  stage  3 
is  retained. 

5.2.  in  the  light  of  the  above  two  observations,  we  seek  a method  of 
reconstruction  which 

uses  the  same  pattern  of  decomposition  as  the  Phase  1 and  2 procedures, 
and  avoids  numerical  difficulties  discussed  earlier, 
is  able  to  reconstruct  at  any  point,  not  just  at  optimality, 
accurately  mirrors  Method  I,  but  gets  around  having  to  save  all 
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(1) 


(2) 

(3) 


extreme  points  generated.  We  accept  it  to  be  reasonable  that  if 

a proposal  B^x^  is  saved,  then  the  corresponding  ^ Is  also  saved. 

When  B . is  throvm  away,  then  the  corresponding  xj  can  be  dispensed 
i— 1 l 

with. 

(4)  is  able  to  utilize  a strategy  akin  to  that  used  in  a standard  simplex 

* 

method,  when  the  optimal  values  of  a set  of  variables  say  x^ , are  known. 


We  suggest  the  following  strategy: 

Method  III: 

* * 

Step  1:  Let  \ and  t m 

c — t m 

* * 

Step  2:  At  stage  t use  At  to  compute  by 

4-i  ■ J k«4-i 

(Recall  (3)  above.) 

If  (t-1)  - 1,  then  STOP. 

Step  3:  Omit  stage  t and  all  subsequent  stages,  put  upper  and  lower 

bounds  on  x , , constraining  these  variables  to  be  close 
— t-i 

* * 
to  , and  solve  the  (t-1)  stage  problem  for  A ^ using 

nested  decomposition. 

Step  4.  t *■  (t-1)  and  GOTO  Step  2. 


The  major  difference  in  Phase  3 between  Method  III  and  Method  II 
is  that  Method  III  uses  the  same  pattern  of  decomposition  as  the  one  used  in 
Phase  2.  The  price  one  pays  for  this  is  the  added  information  that  must  be  saved, 
as  discussed  under  (3)  above  in  this  Section. 
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6.  Conclusions 

We  conclude  by  recounting  some  lessons  we  have  learned  from  this 
effort  for  the  implementation  of  large  scale  LP  algorithms,  and  some  specific 
recommendations  for  nested  decomposition. 

(a)  This  effort  has  benefited  substantially  from  the  numerical  experiences 
of  M.  Aganaglc  [1977]  on  a real  life  model,  and  we  reiterate  the 
Importance  of  gathering  such  practical  experience.  Indeed  this  is  one 
of  the  points  emphasized  in  Dantzig  and  Parikh  [1977],  who  discuss  how 
the  activities  of  algorithms  and  model  development  complement  each  other. 

(b)  at  the  same  time  it  is  clear  that  the  experimental  Information  gathered 
will  be  suspect,  unless  the  experimental  implementation  is  sound  and 

is  designed  to  fail  gracefully,  i.e.  it  does  not  have  to  deal  with  all 
possible  cases,  but  it  must  at  least  provide  reasonably  good  clues 
when  the  algorithm  encounters  difficulties.  This  is  the  theme  of 
Nazareth  [1978a],  where  a set  of  software  aids  for  developing  good 


experimental  implementations  is  described. 

(c)  our  research  demonstrates  that  a fairly  substantial  reformulation  of 
the  nested  decomposition  algorithm  may  be  necessary  before  setting 
about  developing  a robust  user  oriented  version,  which  is  capable  of 
handling  a diverse  set  of  LP  applications. 

(d)  because  error  propagation  may  be  an  inherent  problem  we  have  been 
prompted  to  look  at  other  decompositions,  see  Nazareth  [1978b]. 

(e)  It  is  clear  that  some  extremely  interesting  numerical  problems  arise 
in  nested  decomposition.  Our  results  are  only  a small  contribution, 
and  we  feel  that  a more  complete  error  analysis  and  perturbation  theory  , 
which  we  now  plan  to  undertake,  would  be  a substantial  contribution. 
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There  is  also  a definite  need  for  a sound  experimental  implementation, 

f 

in  order  to  carry  out  further  experimentation.  This  code  should,  for 
example, 

i 

V 

(i)  use  stable  factoring  and  updating  techniques  as  described  in  Reid  [1976]. 

(ii)  use  the  method  of  reconstruction  described  in  Section  5. 

(iii)  use  the  techniques  of  structured  programming,  since  the  data  flow 

in  the  algorithm  is  quite  complex.  Data  flow  and  numeric  consider- 
ations should  be  separated  when  possible. 

(iv)  be  able  to  start  from  an  arbitrary  feasible  (not  necessarily  basic) 
solution,  either  user  supplied  or  developed  from  an  intermediate 
reconstruction. 

(v)  use  a technique  similar  to  that  in  the  GUB  algorithm,  Dantzig  and 
Van  Slyke  [1967]  with  LU  factorization  of  the  basis,  see  Tomlin 
[1974],  to  ensure  that  convexity  rows  are  closely  satisfied. 
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